<Serializable()> Public Class Matrix3DInt16
    Implements IDisposable

    Public data(,,) As Int16
    Public x, y, z As Integer
    Public Sub New()
        ReDim data(0, 0, 0)
        data(0, 0, 0) = 0
        x = 1
        y = 1
        z = 1
    End Sub
    Public Sub New(ByVal x As Integer, ByVal y As Integer, ByVal z As Integer)
        Me.x = x
        Me.y = y
        Me.z = z
        'automatically assigned 0
        ReDim Me.data(0 To x - 1, 0 To y - 1, 0 To z - 1)
    End Sub
    Public Sub New(ByVal data(,,) As Int16)
        Me.x = data.GetUpperBound(0) + 1
        Me.y = data.GetUpperBound(1) + 1
        Me.z = data.GetUpperBound(2) + 1
        ReDim Me.data(0 To x - 1, 0 To y - 1, 0 To z - 1)
        For i As Integer = 0 To x - 1
            For j As Integer = 0 To y - 1
                For k As Integer = 0 To z - 1
                    Me.data(i, j, k) = data(i, j, k)
                Next
            Next
        Next
    End Sub

    Function Minimum() As Int16
        Dim ans As Single = data(0, 0, 0)
        Dim i, j, k As Integer
        For i = 0 To x - 1
            For j = 0 To y - 1
                For k = 0 To z - 1
                    If data(i, j, k) < ans Then
                        ans = data(i, j, k)
                    End If
                Next
            Next
        Next
        If (ans = Maximum()) Then ans = 0
        Return ans
    End Function

    Function Maximum() As Int16
        Dim ans As Single = 0
        Dim i, j, k As Integer
        For i = 0 To x - 1
            For j = 0 To y - 1
                For k = 0 To z - 1
                    If data(i, j, k) > ans Then
                        ans = data(i, j, k)
                    End If
                Next
            Next
        Next
        Return ans
    End Function

    Private Disposed As Boolean = False
    Private Overloads Sub Dispose(ByVal disposing As Boolean)
        ' Check to see if Dispose has already been called.
        If Not Me.Disposed Then
            ' If disposing equals true, dispose all managed 
            ' and unmanaged resources.
            If disposing Then
                ' Dispose managed resources.
                x = Nothing
                y = Nothing
                z = Nothing
                data = Nothing
            End If

            ' Call the appropriate methods to clean up 
            ' unmanaged resources here.
            ' If disposing is false, 
            ' only the following code is executed.
        End If
        Disposed = True
    End Sub
    ' Implement IDisposable.
    ' Do not make this method virtual.
    ' A derived class should not be able to override this method.
    Public Overloads Sub Dispose() Implements IDisposable.Dispose
        Dispose(True)
        ' This object will be cleaned up by the Dispose method.
        ' Therefore, you should call GC.SupressFinalize to
        ' take this object off the finalization queue 
        ' and prevent finalization code for this object
        ' from executing a second time.
        GC.SuppressFinalize(Me)
    End Sub

    Function TransverseSlice(ByVal zslice As Integer) As Int16(,)
        Dim ans(Me.x, Me.y) As Int16
        Dim i, j As Integer
        For i = 0 To Me.x - 1
            For j = 0 To Me.y - 1
                ans(i, j) = Me.data(i, j, zslice)
            Next
        Next
        Return ans
    End Function
    Function SliceToBitmap(ByVal data As Int16(,)) As System.Drawing.Bitmap
        Dim a(,) As Int16 = data
        Dim image As New Drawing.Bitmap(a.GetUpperBound(0) + 1, a.GetUpperBound(1) + 1)
        Dim max As Int16 = Maximum() 'global maximum
        Dim temp As Int16
        For i As Integer = 0 To a.GetUpperBound(0)
            For j As Integer = 0 To a.GetUpperBound(1)
                temp = CShort(a(i, j) * 255 / max)
                image.SetPixel(i, j, System.Drawing.Color.FromArgb(temp, temp, temp))
            Next
        Next
        Return image
    End Function
    Function CoronalSlice(ByVal yslice As Integer) As Int16(,)
        Dim ans(Me.x, Me.z) As Int16
        Dim i, j As Integer
        For i = 0 To Me.x - 1
            For j = 0 To Me.z - 1
                ans(i, Me.z - 1 - j) = Me.data(i, yslice, j)
            Next
        Next
        Return ans
    End Function
    Function SagittalSlice(ByVal xslice As Integer) As Int16(,)
        Dim ans(Me.y, Me.z) As Int16
        Dim i, j As Integer
        For i = 0 To Me.y - 1
            For j = 0 To Me.z - 1
                ans(Me.y - 1 - i, Me.z - 1 - j) = Me.data(xslice, i, j)
            Next
        Next
        Return ans
    End Function
    Public Function Size() As Coordinate
        Return New Coordinate(Me.x, Me.y, Me.z)
    End Function
    '''<Summary>
    '''This function does a 3D Rigid Body Transformation
    '''</Summary>
    '''<param name="ThetaX">Rotation angle in degrees along X axis</param>
    '''<param name="ThetaY">Rotation angle in degrees along Y axis</param>
    '''<param name="ThetaZ">Rotation angle in degrees along Z axis</param>
    '''<param name="TranslateX">Translation along X Axis</param>
    '''<param name="TranslateY">Translation along Y Axis</param>
    '''<param name="TranslateZ">Translation along Z axis</param>
    '''<returns>Transformed 3D Matrix</returns>
    '''<see>http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=7107&objectType=file</see>
    Public Sub RotateAndTranslate(ByVal ThetaX As Single, ByVal ThetaY As Single, ByVal ThetaZ As Single, ByVal TranslateX As Single, ByVal TranslateY As Single, ByVal TranslateZ As Single)
        Dim Xcenter, Ycenter, Zcenter As Single
        Xcenter = Me.x / 2
        Ycenter = Me.y / 2
        Zcenter = Me.z / 2
        ThetaZ = -(ThetaZ * 3.14 / 180)
        ThetaX = -(ThetaX * 3.14 / 180)
        ThetaY = -(ThetaY * 3.14 / 180)

        Dim RotateMatrixX, RotateMatrixY, RotateMatrixZ, FinalRotateMatrix As New twoDMatrixLibrary.Matrix(4, 4)

        RotateMatrixX(0, 0) = 1
        RotateMatrixX(0, 1) = 0
        RotateMatrixX(0, 2) = 0
        RotateMatrixX(0, 3) = 0
        RotateMatrixX(1, 0) = 0
        RotateMatrixX(1, 1) = Math.Cos(ThetaX)
        RotateMatrixX(1, 3) = -1 * Math.Sin(ThetaX)
        RotateMatrixX(1, 3) = 0
        RotateMatrixX(2, 0) = 0
        RotateMatrixX(2, 1) = Math.Sin(ThetaX)
        RotateMatrixX(2, 2) = Math.Cos(ThetaX)
        RotateMatrixX(2, 3) = 0
        RotateMatrixX(3, 0) = 0
        RotateMatrixX(3, 1) = 0
        RotateMatrixX(3, 2) = 0
        RotateMatrixX(3, 3) = 1
        RotateMatrixY(0, 0) = Math.Cos(ThetaY)
        RotateMatrixY(0, 1) = 0
        RotateMatrixY(0, 2) = Math.Sin(ThetaY)
        RotateMatrixY(0, 3) = 0
        RotateMatrixY(1, 0) = 0
        RotateMatrixY(1, 1) = 1
        RotateMatrixY(1, 3) = 0
        RotateMatrixY(1, 3) = 0
        RotateMatrixY(2, 0) = -1 * Math.Sin(ThetaY)
        RotateMatrixY(2, 1) = 0
        RotateMatrixY(2, 2) = Math.Cos(ThetaY)
        RotateMatrixY(2, 3) = 0
        RotateMatrixY(3, 0) = 0
        RotateMatrixY(3, 1) = 0
        RotateMatrixY(3, 2) = 0
        RotateMatrixY(3, 3) = 1
        RotateMatrixZ(0, 0) = Math.Cos(ThetaZ)
        RotateMatrixZ(0, 1) = -1 * Math.Sin(ThetaZ)
        RotateMatrixZ(0, 2) = 0
        RotateMatrixZ(0, 3) = 0
        RotateMatrixZ(1, 0) = Math.Sin(ThetaZ)
        RotateMatrixZ(1, 1) = Math.Cos(ThetaZ)
        RotateMatrixZ(1, 3) = 0
        RotateMatrixZ(1, 3) = 0
        RotateMatrixZ(2, 0) = 0
        RotateMatrixZ(2, 1) = 0
        RotateMatrixZ(2, 2) = 1
        RotateMatrixZ(2, 3) = 0
        RotateMatrixZ(3, 0) = 0
        RotateMatrixZ(3, 1) = 0
        RotateMatrixZ(3, 2) = 0
        RotateMatrixZ(3, 3) = 1

        FinalRotateMatrix = RotateMatrixX * RotateMatrixY * RotateMatrixZ

        Dim TranslationmMatrix As New twoDMatrixLibrary.Matrix(1, 3)
        TranslationmMatrix(0, 0) = TranslateX
        TranslationmMatrix(0, 1) = TranslateY
        TranslationmMatrix(0, 2) = TranslateZ

        Dim temp(,,) As Short
        ReDim temp(x, y, z)
        Dim ii, jj, kk As Integer
        For k As Integer = 0 To z - 1
            For i As Integer = 0 To x - 1
                For j As Integer = 0 To y - 1
                    ii = CInt((FinalRotateMatrix(0, 0) * (i - Xcenter)) + (FinalRotateMatrix(0, 1) * (j - Ycenter)) + (FinalRotateMatrix(0, 2) * (k - Zcenter)) + Xcenter - TranslationmMatrix(0, 0))
                    jj = CInt((FinalRotateMatrix(1, 0) * (i - Xcenter)) + (FinalRotateMatrix(1, 1) * (j - Ycenter)) + (FinalRotateMatrix(1, 2) * (k - Zcenter)) + Ycenter - TranslationmMatrix(0, 1))
                    kk = CInt((FinalRotateMatrix(2, 0) * (i - Xcenter)) + (FinalRotateMatrix(2, 1) * (j - Ycenter)) + (FinalRotateMatrix(2, 2) * (k - Zcenter)) + Zcenter - TranslationmMatrix(0, 2))
                    If ii > 0 And jj > 0 And kk > 0 And ii <= x And jj <= y And kk <= z Then
                        temp(i, j, k) = data(ii - 1, jj - 1, kk - 1)
                    End If
                Next
            Next
        Next
        data = temp
    End Sub
    Public Function RotateAndTranslateFunction(ByVal ThetaX As Single, ByVal ThetaY As Single, ByVal ThetaZ As Single, ByVal TranslateX As Single, ByVal TranslateY As Single, ByVal TranslateZ As Single) As Matrix3DInt16
        Dim tempm As New Matrix3DInt16(data)
        Dim Xcenter, Ycenter, Zcenter As Single
        Xcenter = Me.x / 2
        Ycenter = Me.y / 2
        Zcenter = Me.z / 2
        ThetaZ = -(ThetaZ * 3.14 / 180)
        ThetaX = -(ThetaX * 3.14 / 180)
        ThetaY = -(ThetaY * 3.14 / 180)

        Dim RotateMatrixX, RotateMatrixY, RotateMatrixZ, FinalRotateMatrix As New twoDMatrixLibrary.Matrix(4, 4)


        RotateMatrixX(0, 0) = 1
        RotateMatrixX(0, 1) = 0
        RotateMatrixX(0, 2) = 0
        RotateMatrixX(0, 3) = 0
        RotateMatrixX(1, 0) = 0
        RotateMatrixX(1, 1) = Math.Cos(ThetaX)
        RotateMatrixX(1, 3) = -1 * Math.Sin(ThetaX)
        RotateMatrixX(1, 3) = 0
        RotateMatrixX(2, 0) = 0
        RotateMatrixX(2, 1) = Math.Sin(ThetaX)
        RotateMatrixX(2, 2) = Math.Cos(ThetaX)
        RotateMatrixX(2, 3) = 0
        RotateMatrixX(3, 0) = 0
        RotateMatrixX(3, 1) = 0
        RotateMatrixX(3, 2) = 0
        RotateMatrixX(3, 3) = 1

        RotateMatrixY(0, 0) = Math.Cos(ThetaY)
        RotateMatrixY(0, 1) = 0
        RotateMatrixY(0, 2) = Math.Sin(ThetaY)
        RotateMatrixY(0, 3) = 0
        RotateMatrixY(1, 0) = 0
        RotateMatrixY(1, 1) = 1
        RotateMatrixY(1, 3) = 0
        RotateMatrixY(1, 3) = 0
        RotateMatrixY(2, 0) = -1 * Math.Sin(ThetaY)
        RotateMatrixY(2, 1) = 0
        RotateMatrixY(2, 2) = Math.Cos(ThetaY)
        RotateMatrixY(2, 3) = 0
        RotateMatrixY(3, 0) = 0
        RotateMatrixY(3, 1) = 0
        RotateMatrixY(3, 2) = 0
        RotateMatrixY(3, 3) = 1

        RotateMatrixZ(0, 0) = Math.Cos(ThetaZ)
        RotateMatrixZ(0, 1) = -1 * Math.Sin(ThetaZ)
        RotateMatrixZ(0, 2) = 0
        RotateMatrixZ(0, 3) = 0
        RotateMatrixZ(1, 0) = Math.Sin(ThetaZ)
        RotateMatrixZ(1, 1) = Math.Cos(ThetaZ)
        RotateMatrixZ(1, 3) = 0
        RotateMatrixZ(1, 3) = 0
        RotateMatrixZ(2, 0) = 0
        RotateMatrixZ(2, 1) = 0
        RotateMatrixZ(2, 2) = 1
        RotateMatrixZ(2, 3) = 0
        RotateMatrixZ(3, 0) = 0
        RotateMatrixZ(3, 1) = 0
        RotateMatrixZ(3, 2) = 0
        RotateMatrixZ(3, 3) = 1

        FinalRotateMatrix = RotateMatrixX * RotateMatrixY * RotateMatrixZ

        Dim TranslationmMatrix As New twoDMatrixLibrary.Matrix(1, 3)
        TranslationmMatrix(0, 0) = TranslateX
        TranslationmMatrix(0, 1) = TranslateY
        TranslationmMatrix(0, 2) = TranslateZ

        Dim temp(,,) As Short
        ReDim temp(x, y, z)
        Dim ii, jj, kk As Integer
        For k As Integer = 0 To z - 1
            For i As Integer = 0 To x - 1
                For j As Integer = 0 To y - 1
                    ii = CInt((FinalRotateMatrix(0, 0) * (i - Xcenter)) + (FinalRotateMatrix(0, 1) * (j - Ycenter)) + (FinalRotateMatrix(0, 2) * (k - Zcenter)) + Xcenter - TranslationmMatrix(0, 0))
                    jj = CInt((FinalRotateMatrix(1, 0) * (i - Xcenter)) + (FinalRotateMatrix(1, 1) * (j - Ycenter)) + (FinalRotateMatrix(1, 2) * (k - Zcenter)) + Ycenter - TranslationmMatrix(0, 1))
                    kk = CInt((FinalRotateMatrix(2, 0) * (i - Xcenter)) + (FinalRotateMatrix(2, 1) * (j - Ycenter)) + (FinalRotateMatrix(2, 2) * (k - Zcenter)) + Zcenter - TranslationmMatrix(0, 2))
                    If ii > 0 And jj > 0 And kk > 0 And ii <= x And jj <= y And kk <= z Then
                        temp(i, j, k) = data(ii - 1, jj - 1, kk - 1)
                    End If
                Next
            Next
        Next
        tempm.data = temp
        Return tempm
    End Function
    'Public Function DownSample(ByVal BinSize As Integer)
    '    Dim Min As Integer = CInt(Me.Minimum)
    '    Dim Max As Integer = CInt(Me.Maximum)
    '    Dim NumberOfBins As Integer = (Max - Min) / BinSize + IIf((Max - Min) Mod BinSize = 0, 0, 1)

    '    Dim Temp As New Matrix3DInt16(x, y, z)

    '    For i As Integer = 0 To x - 1
    '        For j As Integer = 0 To y - 1
    '            For k As Integer = 0 To z - 1
    '                If data(i, j, k) > ans Then

    '                End If
    '            Next
    '        Next
    '    Next

    'End Function

    Default Public Property Item(ByVal x As Integer, ByVal y As Integer, ByVal z As Integer) As Int16
        Get
            Return data(x, y, z)
        End Get
        Set(ByVal value As Int16)
            data(x, y, z) = value
        End Set
    End Property
    Default Public Property Item(ByVal xD As Double, ByVal yD As Double, ByVal zD As Double) As Int16
        Get
            Return Me(Convert.ToInt16(Math.Round(xD)), Convert.ToInt16(Math.Round(yD)), Convert.ToInt16(Math.Round(zD)))
        End Get
        Set(ByVal value As Int16)
            Me(Convert.ToInt16(Math.Round(xD)), Convert.ToInt16(Math.Round(yD)), Convert.ToInt16(Math.Round(zD))) = value
        End Set
    End Property
    Public Function ReScale(ByVal Data As Matrix3DInt16, _
    Optional ByVal MinimumValueOfRange As Integer = 0, _
    Optional ByVal MaximumValueOfRange As Integer = 255) As Matrix3DInt16
        Dim ReturnMatrix As New Matrix3DInt16(Data.x, Data.y, Data.z)
        Dim Min As Int16 = Data.Minimum
        Dim Max As Int16 = Data.Maximum
        If Max = Min And Max = 0 Then
            Throw New Exception("Maximum = Minimum while resizing. Please check data and if possible use a higher resolution dataset.")
        End If
        If Min = Max And Max <> 0 Then
            Min = 0
        End If
        For i As Integer = 0 To ReturnMatrix.x - 1
            For j As Integer = 0 To ReturnMatrix.y - 1
                For k As Integer = 0 To ReturnMatrix.z - 1
                    ReturnMatrix(i, j, k) = ((Data(i, j, k) - Min) / (Max - Min)) * (MaximumValueOfRange - MinimumValueOfRange) + MinimumValueOfRange
                Next
            Next
        Next
        Return ReturnMatrix
    End Function

    Sub InsertSlice(ByVal Begin0_End1 As Integer, ByVal C_T_S As String)
        Select Case C_T_S
            Case "C"
                Select Case Begin0_End1
                    Case 0
                        Dim temp As New Matrix3DInt16(Me.x, Me.y + 1, Me.z)
                        For a As Integer = 0 To Me.x - 1
                            For b As Integer = 0 To Me.y - 1
                                For c As Integer = 0 To Me.z - 1
                                    temp(a, b + 1, c) = data(a, b, c)
                                Next
                            Next
                        Next
                        data = temp.data
                    Case 1
                        Dim temp As New Matrix3DInt16(Me.x, Me.y + 1, Me.z)
                        For a As Integer = 0 To Me.x - 1
                            For b As Integer = 0 To Me.y - 1
                                For c As Integer = 0 To Me.z - 1
                                    temp(a, b, c) = data(a, b, c)
                                Next
                            Next
                        Next
                        data = temp.data
                End Select
            Case "T"
                Select Case Begin0_End1
                    Case 0
                        Dim temp As New Matrix3DInt16(Me.x, Me.y, Me.z + 1)
                        For a As Integer = 0 To Me.x - 1
                            For b As Integer = 0 To Me.y - 1
                                For c As Integer = 0 To Me.z - 1
                                    temp(a, b, c + 1) = data(a, b, c)
                                Next
                            Next
                        Next
                        data = temp.data
                    Case 1
                        Dim temp As New Matrix3DInt16(Me.x, Me.y, Me.z + 1)
                        For a As Integer = 0 To Me.x - 1
                            For b As Integer = 0 To Me.y - 1
                                For c As Integer = 0 To Me.z - 1
                                    temp(a, b, c) = data(a, b, c)
                                Next
                            Next
                        Next
                        data = temp.data
                End Select
            Case "S"
                Select Case Begin0_End1
                    Case 0
                        Dim temp As New Matrix3DInt16(Me.x + 1, Me.y, Me.z)
                        For a As Integer = 0 To Me.x - 1
                            For b As Integer = 0 To Me.y - 1
                                For c As Integer = 0 To Me.z - 1
                                    temp(a + 1, b, c) = data(a, b, c)
                                Next
                            Next
                        Next
                        data = temp.data
                    Case 1
                        Dim temp As New Matrix3DInt16(Me.x + 1, Me.y, Me.z)
                        For a As Integer = 0 To Me.x - 1
                            For b As Integer = 0 To Me.y - 1
                                For c As Integer = 0 To Me.z - 1
                                    temp(a, b, c) = data(a, b, c)
                                Next
                            Next
                        Next
                        data = temp.data
                End Select
        End Select
    End Sub

    Public Function MidSliceContrastStretch(Optional ByVal Percentile As Double = 5) As Matrix3DInt16
        Dim Temp As New Matrix3DInt16(x, y, z)
        Dim min As Int16 = Me.Minimum
        Dim max As Int16 = Me.Maximum
        For i As Integer = 0 To x - 1
            For j As Integer = 0 To y - 1
                For k As Integer = 0 To z - 1
                    Temp(i, j, k) = CShort((data(i, j, k) - min) * 255 / (max - min))
                Next
            Next
        Next
        'find histogram
        Dim Histogram(256) As Double
        For i As Integer = 0 To x - 1
            For j As Integer = 0 To y - 1
                'ignore the zero
                If Temp(i, j, Math.Round(z / 2)) <> 0 Then Histogram(Temp(i, j, Math.Round(z / 2))) += 1
            Next
        Next
        'scale histogram
        Dim histmax As Double = 0
        For i As Integer = 0 To 255
            If Histogram(i) > histmax Then
                histmax = Histogram(i)
            End If
        Next
        'scale histogram
        For i As Integer = 0 To 255
            Histogram(i) = Histogram(i) * 100 / histmax
        Next

        'find approx location of 95 and 5 percentiles
        Dim location95 As Integer = 255, location5 As Integer = 0
        For i As Integer = 0 To 255
            If Histogram(i) > Percentile Then
                location5 = i
                Exit For
            End If
        Next
        For i As Integer = 255 To 0 Step -1
            If Histogram(i) > Percentile Then
                location95 = i
                Exit For
            End If
        Next
        location5 = (location5 * (max - min) / 255) + min
        location95 = (location95 * (max - min) / 255) + min
        For i As Integer = 0 To x - 1
            For j As Integer = 0 To y - 1
                For k As Integer = 0 To z - 1
                    Temp(i, j, k) = data(i, j, k)
                    If Temp(i, j, k) > location95 Then Temp(i, j, k) = location95
                    If Temp(i, j, k) < location5 Then Temp(i, j, k) = location5
                Next
            Next
        Next
        Return Temp
    End Function

    Public Function ReScaleToSingle(ByVal Data As Matrix3DInt16, _
        Optional ByVal MinimumValueOfRange As Single = 0, _
        Optional ByVal MaximumValueOfRange As Single = 255) As Matrix3DSingle

        Dim ReturnMatrix As New Matrix3DSingle(Data.x, Data.y, Data.z)
        Dim Min As Single = Data.Minimum
        Dim Max As Single = Data.Maximum

        For i As Integer = 0 To ReturnMatrix.x - 1
            For j As Integer = 0 To ReturnMatrix.y - 1
                For k As Integer = 0 To ReturnMatrix.z - 1
                    If Data(i, j, k) <> 0 Then 'keep the noise same
                        ReturnMatrix(i, j, k) = Convert.ToSingle((Data(i, j, k) - Min) * ((MaximumValueOfRange - MinimumValueOfRange) / (Max - Min)) + MinimumValueOfRange)
                    End If
                Next
            Next
        Next

        Return ReturnMatrix
    End Function

End Class
